Data Ingestion - Intake

Overview

In the last notebook, you learned how to efficiently load data from the Microsoft Planetary Computer platform. If that approach works for you, please proceed to a workflow example. However, there are many other established approaches and platforms (e.g. NASA EarthData, USGS EarthExplorer, Google, and AWS), for loading Landsat data that may not preprocess the data in the exact same way or use a similar API. For this reason, in this notebook we will demonstrate common alternative approaches and techniques for general data access, centered around Intake.

First, to demonstrate simple cases, we will begin by reading in a local files with Pandas and xarray. Then we will demonstrate the combined use of Intake and Dask libraries to efficiently fetch data from a remote server.

Info

A great way to contribute to this cookbook is to create a notebook that focuses on data access from a specific provider.

Prerequisites

Concepts	Importance	Notes
Intro to Landsat	Necessary	Background
Data Ingestion - Planetary Computer	Helpful
Pandas Cookbook	Helpful
xarray Cookbook	Necessary
Intake Quickstart	Helpful
Intake Cookbook	Necessary

Time to learn: 20 minutes

Imports

import pandas as pd
import xarray as xr
import intake
import hvplot.intake
import hvplot.xarray

import warnings
warnings.simplefilter('ignore', FutureWarning) # Ignore warning about the format of epsg codes

Loading local metadata with Pandas

Often, it’s common to have local files that could contain metadata. Therefore, we’ll start with the simple case of loading small local tabular data, such as .csv files into Pandas Dataframes.

Let’s load information about Landsat 5 bands:

landsat5_bands = pd.read_csv('data/landsat5_bands.csv').set_index('Band')
landsat5_bands

	Description	Range (nm)
Band
1	Blue	450-520
2	Green	520-600
3	Red	630-690
4	Near-Infrared	760-900
5	Short Wavelength Infrared 1	1550-1750
6	Thermal Infrared	10400-12500
7	Short Wavelength Infrared 2	2080-2350

…and the relevant subset of Landsat 8 bands:

landsat8_bands = pd.read_csv('data/landsat8_bands.csv').set_index('Band')
landsat8_bands

	Description	Range (nm)
Band
1	Coastal Aerosol	435-451
2	Blue	452-512
3	Green	533-590
4	Red	636-673
5	Near-Infrared	851-879
6	Short Wavelength Infrared 1	1566-1651
7	Short Wavelength Infrared 2	2107-2294

Interestingly, we can see that the range and index of the bands differ between the satellites. For example, ‘Red’ is Band-4 within Landsat 8, but Band-3 within Landsat 5. This is important to keep in mind whenever comparing data across Landsat missions.

Loading data of any size and location using Intake

Even if your current data files are local, earth science data is often complex, remote, and large, so it’s a good idea to become familiar with an approach that can handle data of nearly any complexity, location, and size. For this reason, we will demonstrate a couple different approaches of the Intake library.

First, we will start with a simpler case of reading the same local tabular file as above with Pandas. Then we will progressively cover more complex cases, such as loading multiple files and loading gridded data from a remote location with an Intake catalog.

ds = intake.open_csv('./data/landsat5_bands.csv')

Importantly, the data is not loaded yet - open_csv only accesses the metadata:

ds

csv:
  args:
    urlpath: ./data/landsat5_bands.csv
  description: ''
  driver: intake.source.csv.CSVSource
  metadata: {}

We can also peak at what format the data will be in when it gets loaded using .container:

ds.container

'dataframe'

Now, to get better insight into the data without loading it all in just yet, we can inspect the data using .to_dask(). You could also use .read_chunked() to return an iterator over the chunks (or partitions) of data. Refer to the ‘Reading Data’ Intake documentation for more information about these options.

dd = ds.to_dask()
dd.head()

	Band	Description	Range (nm)
0	1	Blue	450-520
1	2	Green	520-600
2	3	Red	630-690
3	4	Near-Infrared	760-900
4	5	Short Wavelength Infrared 1	1550-1750

dd.info()

<class 'dask.dataframe.core.DataFrame'>
Columns: 3 entries, Band to  Range (nm)
dtypes: object(2), int64(1)

Notice that now running .info() gives us a lot less information than when we loaded the data earlier using pd.read_csv, this is because the Dask representation of the data only loads what it needs, when it is needed for a computation. However, you can always explicitly load all the data into memory using .read() (but be careful to only run this on a dataset that will fit on your computer’s memory). Our sample dataset is very small, so let’s try it:

df = ds.read()
df.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 7 entries, 0 to 6
Data columns (total 3 columns):
 #   Column        Non-Null Count  Dtype 
---  ------        --------------  ----- 
 0   Band          7 non-null      int64 
 1    Description  7 non-null      object
 2    Range (nm)   7 non-null      object
dtypes: int64(1), object(2)
memory usage: 296.0+ bytes

Loading multiple local files

Intake also lets the user load and concatenate data across multiple files in one command

ds = intake.open_csv(['./data/landsat5_bands.csv', './data/landsat8_bands.csv'])

If the data file names allow, this can be more simply expressed using *:

ds = intake.open_csv('./data/landsat*_bands.csv')

df = ds.read()
df.info()

<class 'pandas.core.frame.DataFrame'>
Int64Index: 14 entries, 0 to 6
Data columns (total 3 columns):
 #   Column        Non-Null Count  Dtype 
---  ------        --------------  ----- 
 0   Band          14 non-null     int64 
 1    Description  14 non-null     object
 2    Range (nm)   14 non-null     object
dtypes: int64(1), object(2)
memory usage: 448.0+ bytes

But sometimes there is data encoded in a file name or path that causes concatenated data to lose some important context. In this example, we lose the information about which version of landsat the training was done on (5 vs 8). To keep track of that information, we can use a python format string {variable:type} to specify our path and declare a new field on our data. That field will get populated based on its value in the path.

ds = intake.open_csv('./data/landsat{Landsat:d}_bands.csv')
df = ds.read().set_index(['Landsat', 'Band'])
df

		Description	Range (nm)
Landsat	Band
5	1	Blue	450-520
	2	Green	520-600
	3	Red	630-690
	4	Near-Infrared	760-900
	5	Short Wavelength Infrared 1	1550-1750
	6	Thermal Infrared	10400-12500
	7	Short Wavelength Infrared 2	2080-2350
8	1	Coastal Aerosol	435-451
	2	Blue	452-512
	3	Green	533-590
	4	Red	636-673
	5	Near-Infrared	851-879
	6	Short Wavelength Infrared 1	1566-1651
	7	Short Wavelength Infrared 2	2107-2294

Loading remote data with a catalog

Probably the most versatile and reproducible approach to loading data with Intake is to use a simple Catalog file to declare how the data should be loaded. The catalog lays out how the data should be loaded, defines some metadata, and specifies any patterns in the file path that should be included in the data. Here is an example of an entry in a catalog YAML file:

sources:
  landsat_5_small:
    description: Small version of Landsat 5 Surface Reflectance Level-2 Science Product.
    driver: rasterio
    cache:
      - argkey: urlpath
        regex: 'earth-data/landsat'
        type: file
    args:
      urlpath: 's3://earth-data/landsat/small/LT05_L1TP_042033_19881022_20161001_01_T1_sr_band{band:d}.tif'
      chunks:
        band: 1
        x: 50
        y: 50
      concat_dim: band
      storage_options: {'anon': True}

The urlpath can be a path to a file, list of files, or a path with glob notation. Alternatively, the path can be written as a python style format_string. In the case where the urlpath is a format string (for example, {band:d}), the fields specified in that string will be parsed from the filenames and returned in the data.

Let’s look at what our catalog contains using open_catalog:

cat = intake.open_catalog('./data/catalog.yml')
list(cat)

['landsat_5_small',
 'landsat_8_small',
 'landsat_5',
 'landsat_8',
 'google_landsat_band',
 'amazon_landsat_band',
 'fluxnet_daily',
 'fluxnet_metadata',
 'seattle_lidar']

Here, we see a list of datasets specified by the catalog.yml file. Lets’s look at a description and the urlpath for one of them:

landsat_5 = cat.landsat_5_small

landsat_5.description

'Small version of Landsat 5 Surface Reflectance Level-2 Science Product.'

landsat_5.container

'xarray'

From .container we can see that ultimately the data will be an xarray object

landsat_5.urlpath

's3://earth-data/landsat/small/LT05_L1TP_042033_19881022_20161001_01_T1_sr_band*.tif'

From the urlpath we can see that this dataset is hosted remotely on Amazon S3

Let’s preview this dataset using .to_dask().

landsat_5 = cat.landsat_5_small
landsat_5.to_dask()

2023-04-28 07:47:13,197 - intake - WARNING - cache.py:_download:L261 - Cache progress bar in a notebook requires ipywidgets to be installed: conda/pip install ipywidgets

2023-04-28 07:47:13,199 - intake - WARNING - cache.py:_download:L261 - Cache progress bar in a notebook requires ipywidgets to be installed: conda/pip install ipywidgets

2023-04-28 07:47:13,528 - intake - WARNING - cache.py:_download:L261 - Cache progress bar in a notebook requires ipywidgets to be installed: conda/pip install ipywidgets

2023-04-28 07:47:13,617 - intake - WARNING - cache.py:_download:L261 - Cache progress bar in a notebook requires ipywidgets to be installed: conda/pip install ipywidgets

2023-04-28 07:47:13,688 - intake - WARNING - cache.py:_download:L261 - Cache progress bar in a notebook requires ipywidgets to be installed: conda/pip install ipywidgets

2023-04-28 07:47:13,760 - intake - WARNING - cache.py:_download:L261 - Cache progress bar in a notebook requires ipywidgets to be installed: conda/pip install ipywidgets

<xarray.DataArray (band: 6, y: 300, x: 300)>
dask.array<concatenate, shape=(6, 300, 300), dtype=float64, chunksize=(1, 50, 50), chunktype=numpy.ndarray>
Coordinates:
  * band     (band) int64 1 2 3 4 5 7
  * y        (y) float64 4.309e+06 4.309e+06 4.309e+06 ... 4.264e+06 4.264e+06
  * x        (x) float64 3.324e+05 3.326e+05 3.327e+05 ... 3.771e+05 3.772e+05
Attributes:
    transform:      (150.0, 0.0, 332325.0, 0.0, -150.0, 4309275.0)
    crs:            +init=epsg:32611
    res:            (150.0, 150.0)
    is_tiled:       0
    nodatavals:     (nan,)
    scales:         (1.0,)
    offsets:        (0.0,)
    AREA_OR_POINT:  Area

From this information we can see that there are three dimensions with coordinates x, y, and band. This is satellite data where the x and y are spatial earth coordinates and the band refers to the wavelength range of light used to capture the satellite image.

Visualizing the data

Often, data ingestion involves quickly visualizing your raw data to get a sense that things are proceeding accordingly. You have already seen how we can easily use hvPlot for bare xarray data, but hvPlot also provides interactive plotting commands for Intake, Pandas, Dask, and GeoPandas objects. We’ll look more closely at hvPlot and its options in later tutorials.

For now, let’s just plot out our different band images as columns (using col='band'). We can also set rasterize=True to use Datashader (another HoloViz tool) to render large data into a 2D histogram, where every array cell counts the data points falling into that pixel, as set by the resolution of your screen.

landsat_5.hvplot.image(col='band', cmap='viridis', xaxis=False, yaxis=False, colorbar=False, rasterize=True)

Hint: Try zooming in on one plot to see that they are all linked!

If we increase the complexity a bit and consider the likely situation where you are constantly loading and visualizing a similar type of data using an Intake catalog, you can instead just declare plotting parameters directly into the catalog!

Here is the relevant part of our catalog.yml:

metadata:
  plots:
    band_image:
      kind: 'image'
      x: 'x'
      y: 'y'
      groupby: 'band'
      rasterize: True
      width: 400
      dynamic: False

landsat_5.hvplot.band_image()

Excellent! Since we didn’t specify that the images should to be plotted as columns, hvPlot conveniently gives us a widget to scrub through them.

Side-note: Drilling down for data access

As part of accessing data, it’s important to understand how to drill down to different representations of our data that certain processing steps might require.

Numpy Array

Earlier, we saw that if the selected data are small enough, we can load return all the data into one in-memory xarray DataArray container. Once in an DataArray object, the data can be easily manipulated and visualized. Under the hood, this xarray object is essentially comprised of NumPy arrays with the addition of coordinate labels (and other superpowers). However, sometimes machine learning tools (e.g. scikit-learn) only accept NumPy arrays as input. These arrays are accessible to us from the xarray DataArray objects by using the .values attribute.

npa = cat.landsat_5_small.read().values
type(npa)

numpy.ndarray

npa.shape

(6, 300, 300)

npa[:3, :3, :3]

array([[[ 640.,  842.,  864.],
        [ 796.,  774.,  707.],
        [ 975.,  707.,  908.]],

       [[ 810., 1096., 1191.],
        [1096., 1048.,  905.],
        [1286., 1000., 1286.]],

       [[1007., 1345., 1471.],
        [1260., 1259., 1175.],
        [1555., 1175., 1555.]]])

Summary

Before accessing any data, it’s a good idea to start by learning about the context and details of the dataset. This will give you the intuition to make informed decisions as you form a processing and analysis pipeline. Once you have identified a dataset to use, tailor the data access approach given the features of the data, such as size and location. If your data is local and small, a reliable approach is to access the data with either Pandas or xarray, depending on whether the data is tabular or gridded. As earth science data is often large, multidimensional, and on a remote server, a powerful and increasingly common approach is to use the combination of Intake, Dask, and xarray. Once you have accessed data, visualize it with hvPlot to ensure that it matches your expectations.

What’s next?

After loading the data, you may need to process it in some way before it can be fed into your analysis workflow. We will cover an example of such processing in the next section.

Resources and references

Authored/adapted by Demetris Roumis circa Dec, 2022
The Landsat 8 banner image is from NASA